# ======================================================================
# Copyright CERFACS (March 2019)
# Contributor: Adrien Suau (adrien.suau@cerfacs.fr)
#
# This software is governed by the CeCILL-B license under French law and
# abiding  by the  rules of  distribution of free software. You can use,
# modify  and/or  redistribute  the  software  under  the  terms  of the
# CeCILL-B license as circulated by CEA, CNRS and INRIA at the following
# URL "http://www.cecill.info".
#
# As a counterpart to the access to  the source code and rights to copy,
# modify and  redistribute granted  by the  license, users  are provided
# only with a limited warranty and  the software's author, the holder of
# the economic rights,  and the  successive licensors  have only limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using, modifying and/or  developing or reproducing  the
# software by the user in light of its specific status of free software,
# that  may mean  that it  is complicated  to manipulate,  and that also
# therefore  means that  it is reserved for  developers and  experienced
# professionals having in-depth  computer knowledge. Users are therefore
# encouraged  to load and  test  the software's  suitability as  regards
# their  requirements  in  conditions  enabling  the  security  of their
# systems  and/or  data to be  ensured and,  more generally,  to use and
# operate it in the same conditions as regards security.
#
# The fact that you  are presently reading this  means that you have had
# knowledge of the CeCILL-B license and that you accept its terms.
# ======================================================================

r"""Oracles needed to solve the wave equation in 1D with Neumann boundary conditions.

Module content
==============

This module implements the oracles needed to solve the 1-dimensional wave equation.

.. warning::

   For the moment the oracles are split in 3: there is an oracle M for the
   "permutation", an oracle V for the amplitude and a third oracle S for the sign.
   These three oracles will soon be reunited into one big oracle, as described in
   every paper on the subject. This regrouping will allow us to re-use some
   computations that are common to some of the oracles.

Each oracle is implemented as a function taking no argument and returning a QRoutine.
The QRoutines construction has been encapsulated into functions in order to avoid
import-time computations.


Oracle construction explanations
================================

Graphs associated with the chosen discretisation and boundary conditions
------------------------------------------------------------------------

This section explains how the oracles have been constructed.

Following :cite:`1711.05394v1`, we should begin by discretising the 1D line and creating
the associated graph. Let :math:`N_d` be the number of discretisation points used. In
:cite:`1711.05394v1`, the authors ignored the two boundary points in the discretisation
as their value is constant through time. Because of this, we will mostly use
:math:`N_c = N_d - 2`, the number of *considered* points.

The graph associated with the discretisation and the Neumann boundary condition is
defined as :math:`G_N = (V, E_N)` such that:

- :math:`V = \left\{ 0, 1, \dots{}, N_c-1\right\}`
- :math:`E_N = \left\{ (i,j) \middle| (i,j)\in V^2 \text{ and } i = j \pm 1 \right\}`.

  In other words, an edge :math:`(i,j)` is in :math:`E_N` if and only if the vertices
  :math:`i` and :math:`j` are neighbours.

The graphical representation of the graph :math:`G_N` for :math:`N_d` discretisation
points is drawn below.

.. admonition:: Graphical representation of :math:`G_N` for the 1-dimensional case \
   with Neumann boundary condition:

   .. image:: /images/non_oriented_neumann_graph.png

.. admonition:: Note on Dirichlet boundary condition

   The graph :math:`G_D = (V, E_D)` associated with the discretisation and **Dirichlet**
   boundary condition is similar to :math:`G_N` except that the 2 vertices at the
   boundary have self-loops. Mathematically,

   .. math::

      E_D = \left\{ (i,j) \middle| (i,j)\in V^2 \text{ and } i = j \pm 1 \right\} \cup
            \left\{ (0,0), (N_c-1, N_c-1) \right\}

   In other words, an edge :math:`(i,j)` is in :math:`E_D` if and only if the vertices
   :math:`i` and :math:`j` are neighbours or :math:`i = j = 0` or :math:`i = j = N_c-1`.

   .. image:: /images/non_oriented_dirichlet_graph.png

From graphs to matrices
-----------------------

The non-oriented graphs defined in the previous sections are used to define the matrix
:math:`B`, used in turn to define the matrix :math:`H`, i.e. the Hamiltonian matrix we
will have to encode in oracles.

In order to construct the matrix :math:`B`, the graph :math:`G_N` (resp. :math:`G_D`)
should be oriented and each vertex and edge should be numbered. Depending on the
chosen orientation and numbering, the expression of :math:`B` and :math:`H` will change
but :math:`B^\dagger B` will not. As explained in :cite:`1711.05394v1`, the quantity
that we will end up simulating is :math:`B^\dagger B` [#bdagb]_, and so the orientation
and numbering of :math:`G_N` (resp. :math:`G_D`) is not constrained.

Because our final goal is to encode the matrix :math:`H` (created from :math:`B`, which
is in turn created from :math:`G_N`), we will choose a orientation and a numbering that
result in a matrix :math:`H` that can be concisely described [#consicelydescribed]_.

.. admonition:: Numbering and ordering for Neumann boundary conditions

   .. image:: /images/oriented_neumann_graph_numbered.png

.. admonition:: Numbering and ordering for Dirichlet boundary conditions

   .. image:: /images/oriented_dirichlet_graph_numbered.png

The matrix :math:`B`
^^^^^^^^^^^^^^^^^^^^

The matrix :math:`B`, of size
:math:`\vert V \vert \times \vert E \vert \equiv N_c \times (N_c+1) \
\equiv (N_d-2) \times (N_d-1)`, is defined by:

.. math::

   B_{ij} =
   \left\{
   \begin{array}{lr}
     \sqrt{W_j} & \text{if } j \text{ is a self-loop of }i \\
     \sqrt{W_j} & \text{if } j \text{ is an edge with } i \text{ as source} \\
    -\sqrt{W_j} & \text{if } j \text{ is an edge with } i \text{ as sink} \\
     0          & \text{otherwise}
   \end{array}
   \right.

In the definition above, the exact role of :math:`W_j` is still unclear. For the moment
it is enough to set :math:`W_j = 1`.

For Dirichlet boundary condition, the matrix :math:`B` is

.. math::

   \newcommand{\Bmatrix}{%
       \begin{matrix}%
       1      & 1      & 0      & \cdots & 0      \\%
       0      & -1     & 1      & \ddots & \vdots \\%
       \vdots & \ddots & \ddots & \ddots & 0      \\%
       0      & \cdots & 0      & -1     & 1      \\%
       \end{matrix}%
   }
   B = \left.\left(
   \vphantom{\Bmatrix}
   \smash{\overbrace{\Bmatrix}^{N_c+1}}\right) \right\} N_c


The matrix :math:`H`
^^^^^^^^^^^^^^^^^^^^

The Hamiltonian :math:`H` we need to simulate to solve the wave equation is given by

.. math::

   H = \begin{pmatrix} 0 & B & 0 \\ B^\dagger & 0 & 0 \\ 0 & 0 & 0 \\ \end{pmatrix}.

The last bottom-lines and right-columns of :math:`B` are full of :math:`0` and are just
here to fill :math:`H` until it is a :math:`\left( 2^q, 2^q \right)` matrix.

Clearly viewing the structure of :math:`H` will be necessary for the oracle construction
part. In order to ease the understanding, the full structure of :math:`H` is recalled
below:

.. math::

   \newcommand{\Hmatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & 0      & 1      & 1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        & \vdots & 0      & -1     & 1      & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      & 0      & \cdots & 0      & -1     & 1      & \vdots &        &        &        &        & \vdots \\%
           1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \vdots &        &        &        &        & \vdots \\%
           1      & -1     & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & 1      & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots & \ddots & \ddots & -1     & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & 0      & 1      & 0      & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0        %
       \end{matrix}%
   }%
   \newcommand{\FirstZeroMatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & 0      \\%
           \vdots &        &        & \vdots \\%
           \vdots &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
   \newcommand{\BTmatrix}{%
       \begin{matrix}%
           1      & 0      & \cdots & 0      \\%
           1      & -1     & \ddots & \vdots \\%
           0      & 1      & \ddots & 0      \\%
           \vdots & \ddots & \ddots & -1     \\%
           0      & \cdots & 0      & 1      \\%
       \end{matrix}%
   }%
   \newcommand{\SecondZeroMatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        & \vdots \\%
           \vdots &        &        &        & \vdots \\%
           \vdots &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
   \newcommand{\RightTopZeroMatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
       \end{matrix}%
   }%
   \newcommand{\RightMidZeroMatrix}{%
       \begin{matrix}%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
   \newcommand{\BotLeftZeroMatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & \cdots \\%
           \vdots &        &        &        \\%
           \vdots &        &        &        \\%
           \vdots &        &        &        \\%
           \vdots &        &        &        \\%
           0      & \cdots & \cdots & \cdots \\%
       \end{matrix}%
   }%
   \newcommand{\BotMidZeroMatrix}{%
       \begin{matrix}%
           \cdots & \cdots & \cdots & \cdots & 0      \\%
                  &        &        &        & \vdots \\%
                  &        &        &        & \vdots \\%
                  &        &        &        & \vdots \\%
                  &        &        &        & \vdots \\%
           \cdots & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
   \newcommand{\BotRightZeroMatrix}{%
       \begin{matrix}%
           0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
    H = \begin{array}{r}
           \begin{array}{cccp{5pt}}
             \overbrace{\hphantom{\BTmatrix}}^{N_c}\!\! &
             \overbrace{\hphantom{\Bmatrix}}^{N_c+1}\!\! &
             \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q - \left( 2N_c + 1 \right)}&
           \end{array}\\[0pt]
           \begin{array}{r}
             N_c                          \left\{\vphantom{\Bmatrix}\right.      \\[5pt]
             N_c + 1                      \left\{\vphantom{\BTmatrix}\right.     \\[5pt]
             2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
           \end{array}
           \left( \Hmatrix \right)
         \end{array}




.. [#bdagb] See :cite:`1711.05394v1` for the full explanation.
.. [#consicelydescribed] This means that :math:`H` has a well defined structure which
   can be described with a few mathematical operations.

From matrices to oracles
------------------------

Presentation of the different oracles needed
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The simulation algorithm used needs three oracles (that can be either given separately
or as one big oracle):

#. :math:`M` is the "permutation" oracle. It takes as input a quantum state
   :math:`\vert x \rangle` and returns in another quantum state the column index of
   the first non-zero element in the row :math:`x`.

   .. math::

      M \vert x \rangle \vert 0 \rangle \mapsto \vert x \rangle \vert m(x) \rangle

   .. tip::
      The :math:`M` oracle does not need to encode a :math:`1`\ -sparse permutation
      matrix. If some rows are filled with zeros (i.e. no non-zero element) then the
      :math:`M` oracle can return any value. This property is used to reduce the
      number of gates needed in some oracles.

#. :math:`V` is the "amplitude" oracle. It takes as input a quantum state
   :math:`\vert x \rangle` and returns in another quantum state the absolute value of
   the first non-zero element in the row :math:`x`: :math:`\vert w(x) \vert`.

   .. math::

      V \vert x \rangle \vert 0 \rangle \mapsto \vert x \rangle \vert \vert w(x)\vert
      \rangle

   .. note::
      If the given row index :math:`\vert x \rangle` has no non-zero element, the oracle
      :math:`V` should return :math:`\vert 0\rangle`.

#. :math:`S` is the "sign" oracle. It takes as input a quantum state
   :math:`\vert x \rangle` and returns in another quantum state the sign
   of the first non-zero element in the row :math:`x`: :math:`s(w(x)) = \left\{  \
   \begin{array}{lr} 0 & w(x) \geqslant 0 \\ 1 & w(x) < 0 \end{array} \
   \right.`.

   .. math::

      S \vert x \rangle \vert 0 \rangle \mapsto \vert x \rangle \vert s(w(x)) \rangle


Decomposition of :math:`H`
^^^^^^^^^^^^^^^^^^^^^^^^^^

The algorithms used to simulate :math:`H` need to have access to a decomposition of
:math:`H` in the form of a sum of :math:`1`\ -sparse Hamiltonian matrices.

Let

.. math::

   H = H_1 + H_2

with:

.. math::

   \newcommand{\HmatrixOne}{%
       \begin{matrix}%
           0      & \cdots & \cdots & 0      & 1      & 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        & \vdots & 0      & -1     & \ddots &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      & 0      & \cdots & 0      & -1     & 0      & \vdots &        &        &        &        & \vdots \\%
           1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \vdots &        &        &        &        & \vdots \\%
           0      & -1     & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        & \ddots & -1     & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
    H_1 = \begin{array}{r}
           \begin{array}{cccp{5pt}}
             \overbrace{\hphantom{\BTmatrix}}^{N_c}\!\! &
             \overbrace{\hphantom{\Bmatrix}}^{N_c+1}\!\! &
             \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q - \left( 2N_c + 1 \right)}&
           \end{array}\\[0pt]
           \begin{array}{r}
             N_c                          \left\{\vphantom{\Bmatrix}\right.      \\[5pt]
             N_c + 1                      \left\{\vphantom{\BTmatrix}\right.     \\[5pt]
             2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
           \end{array}
           \left( \HmatrixOne \right)
         \end{array}

.. math::

   \newcommand{\HmatrixTwo}{%
       \begin{matrix}%
           0      & \cdots & \cdots & 0      & 0      & 1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        & \vdots & \vdots &        & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & 0      & 1      & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \vdots &        &        &        &        & \vdots \\%
           1      & \ddots &        & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \ddots & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & 0      & 1      & 0      & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           \vdots &        &        &        &        &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
           0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
       \end{matrix}%
   }%
    H_2 = \begin{array}{r}
           \begin{array}{cccp{5pt}}
             \overbrace{\hphantom{\BTmatrix}}^{N_c}\!\! &
             \overbrace{\hphantom{\Bmatrix}}^{N_c+1}\!\! &
             \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q - \left( 2N_c + 1 \right)}&
           \end{array}\\[0pt]
           \begin{array}{r}
             N_c                          \left\{\vphantom{\Bmatrix}\right.      \\[5pt]
             N_c + 1                      \left\{\vphantom{\BTmatrix}\right.     \\[5pt]
             2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
           \end{array}
           \left( \HmatrixTwo \right)
         \end{array}


Oracle implementation
^^^^^^^^^^^^^^^^^^^^^

Implementation details are discussed in the documentation of each function implementing
an oracle.



List of possible optimisations that are not implemented at the moment:
----------------------------------------------------------------------

#. The compare function is optimised to apply gates only when needed (i.e. only when
   a bit of the given input is set to 1). In several oracles, the values we compare to
   are not fixed and can be chosen so that they reduce the number of gates.

   .. admonition:: Example

      For the function :func:`~.permutation_oracle_2`, the compare could be performed
      with `discretisation_points_number - 1` **if** the subtraction is guaranteed to
      behave well with underflow.

"""
from qat.lang.AQASM.gates import X, CNOT, CCNOT
from qat.lang.AQASM.routines import QRoutine
from qat.lang.AQASM.misc import build_gate
from qat.lang.AQASM.gates import AbstractGate

from qaths.routines.logical import logical_or


@build_gate(
    "oracle_dirichlet1_1d_wave_equation", [int, int], arity=lambda n, discr: 2 * n + 2
)
def get_oracle_dirichlet1_1d_wave_equation(
    n: int, discretisation_points_number: int
) -> QRoutine:
    r"""Return the oracle for the first matrix in the decomposition.

    The returned routine needs :math:`2n + 2` qubits organised as follow: ::

        |        |x>        |        |m>        |  |v>  |  |s>  |
        |   .   .   .   .   |   .   .   .   .   |   .   |   .   |
        |   0          n-1  |   n         2*n-1 |  2*n  | 2*n+1 |

    The quantum registers :math:`\ket{m}`, :math:`\ket{v}` and :math:`\ket{s}` should
    all be initialised to :math:`\ket{0}`.

    This oracle is split in 3 different oracles: an oracle to encode the positions of
    the non-zero entries (called the permutation oracle), a second oracle to encode
    the absolute values of the weights of the non-zero entries (called the weights
    oracle) and a last oracle encoding the signs of the non-zero entries (the sign
    oracle).

    The returned QRoutine expects 3 AbstractGates:
    - "add_const"
    - "sub_const"
    - "compare_const"

    The permutation oracle encodes the permutation

    .. math::
       \newcommand{\mbf}[1]{\mathbf{#1}}
       \newcommand{\HmatrixOnePerm}{%
           \begin{matrix}%
               0      & \cdots & \cdots & 0      & 1      & 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               \vdots &        &        & \vdots & 0      & \ddots & \ddots &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
               0      & \cdots & \cdots & 0      & 0      & \cdots & 0      & 1      & 0      & \vdots &        &        &        &        & \vdots \\%
               1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \vdots &        &        &        &        & \vdots \\%
               0      & \ddots & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        & \ddots & 1      & 0      &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               0      & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               0      & \cdots & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               \vdots &        &        &        &        & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        &        & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        &        &        & \ddots & \mbf{1}& 0      &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        &        &        &        & 0      & \mbf{1}& \ddots &        &        &        & \vdots \\%
               0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \mbf{1}& 0      & 0      & \cdots & 0      \\%
           \end{matrix}%
       }%
       H_1^M = \begin{array}{r}
               \begin{array}{cccp{5pt}}
                 \overbrace{\hphantom{\BTmatrix}}^{N_c}&
                 \overbrace{\hphantom{\Bmatrix}}^{N_c+1}&
                 \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q -
                     \left( 2N_c + 1 \right)} &
               \end{array}\\[0pt]
               \begin{array}{r}
                 N_c                          \left\{\vphantom{\Bmatrix}\right.  \\[5pt]
                 N_c + 1                      \left\{\vphantom{\BTmatrix}\right. \\[5pt]
                 2^q-\left(2N_c + 1\right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
               \end{array}
               \left( \HmatrixOnePerm \right)
             \end{array}

    by applying the operation

    .. math::

       M_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
       \left\{
       \begin{array}{lr}
       \vert x \rangle \vert x + N_c \rangle & \text{if } x < N_c \\
       \vert x \rangle \vert x - N_c \rangle & \text{else} \\
       \end{array}.
       \right.



    The weights oracle encodes the weights given by the following vector:

    .. math::

       \newcommand{\mbf}[1]{\mathbf{#1}}
       V_1 = \begin{array}{r}
               \begin{array}{r}
                 N_c                          \left\{\vphantom{\begin{array}{l}1
                   \\\vdots \\1\\ \end{array}}\right. \\
                 N_c + 1                      \left\{\vphantom{\begin{array}{l}1
                   \\\vdots \\1\\ \end{array}}\right. \\
                 2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}0
                   \\\vdots \\0\\ \end{array}}\right. \\
               \end{array}
               \begin{pmatrix}
               1 \\
               \vdots \\
               1 \\
               1 \\
               \vdots \\
               1 \\
               0 \\
               \vdots \\
               0
               \end{pmatrix}
             \end{array}

    by applying the operation

    .. math::

       V_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
       \left\{
       \begin{array}{lr}
       \vert x \rangle \vert 1 \rangle & \text{if } x < 2N_c \\
       \vert x \rangle \vert 0 \rangle & \text{else} \\
       \end{array}.
       \right.



    The signs oracle encodes the signs given by the following vector:

    .. math::

       \newcommand{\mbf}[1]{\mathbf{#1}}
       S_1 = \begin{array}{r}
               \begin{array}{r}
                 N_c                          \left\{\vphantom{\begin{array}{l}+\\-
                   \\\vdots \\-\\ \end{array}}\right. \\
                 N_c + 1                      \left\{\vphantom{\begin{array}{l}+\\-
                   \\\vdots \\-\\ \end{array}}\right. \\
                 2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}+/-
                   \\\vdots \\+/-\\ \end{array}}\right.\\
               \end{array}
               \begin{pmatrix}
               + \\
               - \\
               \vdots \\
               - \\
               + \\
               - \\
               \vdots \\
               - \\
               +/- \\
               \vdots \\
               +/-
               \end{pmatrix}
             \end{array}

    by applying the operation

    .. math::

       S_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
       \left\{
       \begin{array}{lr}
       \vert x \rangle \vert 0 \rangle & \text{if } (x = 0) \lor (x = N_c)\\
       \vert x \rangle \vert 1 \rangle & \text{else} \\
       \end{array}.
       \right.

    .. note::

       The :math:`+` sign is encoded as :math:`\vert 0 \rangle`, the :math:`-` sign
       as :math:`\vert 1 \rangle` and :math:`+/-` can be any of :math:`+` or :math:`-`.


    :param n: Size of the :math:`\ket{x}` and :math:`\ket{m}` quantum
        registers.
    :param discretisation_points_number: Number of discretisation points used to
        discretise the wave equation.
    :return: a QRoutine implementing the oracle of the first matrix of the
        decomposition.
    """
    rout = QRoutine()
    x = rout.new_wires(n)
    m = rout.new_wires(n)
    v = rout.new_wires(1)
    s = rout.new_wires(1)
    a = rout.new_wires(4)
    rout.set_ancillae(a)

    # Declaring an AbstractGate that will compare the value stored in a quantum
    # register to a constant value.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to compare.
    # comparator = AbstractGate(
    #     "compare_const", [int, int], arity=lambda q_size, rhs: 2 * q_size
    # )
    comparator = AbstractGate(
        "compare_const", [int, int], arity=lambda q_size, rhs: q_size + 1
    )

    # Declaring an AbstractGate that will add a constant value to the value stored in a
    # quantum register.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to add.
    adder = AbstractGate("add_const", [int, int], arity=lambda q_size, rhs: q_size)

    # Declaring an AbstractGate that will subtract a constant value to the value stored
    # in a quantum register.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to subtract.
    subtractor = AbstractGate("sub_const", [int, int], arity=lambda q_size, rhs: q_size)

    #########################
    # 1. Permutation oracle #
    #########################
    # Copy the value stored in |x> to |m>.
    for i in range(n):
        rout.apply(CNOT, x[i], m[i])

    # Compare |x> with a generation-time value.
    rout.protected_apply(comparator(n, discretisation_points_number - 2), x, a[0])
    # Perform the addition when (|x> < discretisation_points_number-2)
    rout.apply(adder(n, discretisation_points_number - 2).ctrl(), a[0], m)
    # Else, perform a subtraction as the inverse of an addition.
    rout.protected_apply(X, a[0])
    rout.apply(subtractor(n, discretisation_points_number - 2).ctrl(), a[0], m)
    rout.protected_apply(X, a[0])

    # /!\ DO NOT uncompute the comparison as we can re-use it for the sign oracle part.
    # /!\ The compare operation will be uncomputed after the sign oracle.
    # rout.protected_apply(
    #     comparator(n, discretisation_points_number - 2).dag(), x, a[0], m[:-1]
    # )

    #####################
    # 2. Weights oracle #
    #####################
    rout.apply(comparator(n, 2 * (discretisation_points_number - 2)), x, v)

    ###################
    # 3. Signs oracle #
    ###################
    # |a[1]> := (|x>  <  discretisation_points_number-1)
    rout.protected_apply(comparator(n, discretisation_points_number - 1), x, a[1])

    # |a[0]> := ! (|x>  <  discretisation_points_number-2)
    # /!\ the line below has already been computed by the permutation oracle /!\
    # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    # rout.protected_apply(
    #     comparator(n, discretisation_points_number - 2), x, a[0], m[:-1]
    # )
    # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    # Still, we need the negated value.
    rout.protected_apply(X, a[0])

    # |a[2]> := |a[0]> && |a[1]>  <==>  (|x> == discretisation_points_number-2)
    rout.protected_apply(CCNOT, a[0], a[1], a[2])
    # |a[3]> := (|x>  <  1)  <==>  (|x> == 0)
    rout.protected_apply(comparator(n, 1), x, a[3])

    # |s>    := ! (|a[2]> || |a[3]>)
    rout.apply(logical_or(), a[2], a[3], s)
    rout.apply(X, s)

    # Uncompute everything except |s>
    rout.protected_apply(comparator(n, 1).dag(), x, a[3])
    rout.protected_apply(CCNOT.dag(), a[0], a[1], a[2])
    rout.protected_apply(X.dag(), a[0])
    rout.protected_apply(comparator(n, discretisation_points_number - 1).dag(), x, a[1])
    # /!\ Do not forget the uncomputation below /!\
    rout.protected_apply(comparator(n, discretisation_points_number - 2).dag(), x, a[0])

    return rout


@build_gate(
    "oracle_dirichlet2_1d_wave_equation", [int, int], arity=lambda n, discr: 2 * n + 2
)
def get_oracle_dirichlet2_1d_wave_equation(
    n: int, discretisation_points_number: int
) -> QRoutine:
    r"""Return the oracle for the second matrix in the decomposition.

    The returned routine needs :math:`2n + 2` qubits organised as follow: ::

        |        |x>        |        |m>        |  |v>  |  |s>  |
        |   .   .   .   .   |   .   .   .   .   |   .   |   .   |
        |   0          n-1  |   n         2*n-1 |  2*n  | 2*n+1 |

    The quantum registers :math:`\ket{m}`, :math:`\ket{v}` and :math:`\ket{s}` should
    all be initialised to :math:`\ket{0}`.

    This oracle is split in 3 different oracles: an oracle to encode the positions of
    the non-zero entries (called the permutation oracle), a second oracle to encode
    the absolute values of the weights of the non-zero entries (called the weights
    oracle) and a last oracle encoding the signs of the non-zero entries (the sign
    oracle).


    The permutation oracle encodes the permutation

    .. math::

       \newcommand{\mbf}[1]{\mathbf{#1}}
       \newcommand{\HmatrixTwo}{%
           \begin{matrix}%
               0      & \cdots & \cdots & 0      & 0      & 1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        & \vdots & \vdots &        & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
               0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & 0      & 1      & 0      &        &        &        &        & \vdots \\%
               0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \mbf{1}&        &        &        &        & \vdots \\%
               1      & \ddots &        & \vdots & \vdots &        &        &        & \vdots & 0      &        &        &        &        & \vdots \\%
               0      & \ddots & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               0      & \cdots & 0      & 1      & 0      & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               0      & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
               \vdots &        &        &        & \ddots & \ddots & \ddots &        & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        &        & \ddots & \ddots & 0      & 0      &        &        &        &        & \vdots \\%
               \vdots &        &        &        &        &        &        & \ddots & \mbf{1}& 0      & \ddots &        &        &        & \vdots \\%
               0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & \mbf{1}& 0      & 0      & \cdots & \cdots & 0      \\%
           \end{matrix}%
       }%
        H_2^M = \begin{array}{r}
               \begin{array}{cccp{5pt}}
                 \overbrace{\hphantom{\BTmatrix}}^{N_c} &
                 \overbrace{\hphantom{\Bmatrix}}^{N_c+1}&
                 \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q-\left(2N_c + 1\right)}&
               \end{array}\\[0pt]
               \begin{array}{r}
                 N_c                          \left\{\vphantom{\Bmatrix}\right.  \\[5pt]
                 N_c + 1                      \left\{\vphantom{\BTmatrix}\right. \\[5pt]
                 2^q - \left(2N_c+1\right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
               \end{array}
               \left( \HmatrixTwo \right)
             \end{array}

    by applying the operation

    .. math::

       M_{H_2}\vert x \rangle \vert 0 \rangle \mapsto
       \left\{
       \begin{array}{lr}
       \vert x \rangle \vert x + (N_c + 1) \rangle & \text{if } x < (N_c + 1) \\
       \vert x \rangle \vert x - (N_c + 1) \rangle & \text{else} \\
       \end{array}.
       \right.



    The weights oracle encodes the weights given by the following vector:

    .. math::

       \newcommand{\mbf}[1]{\mathbf{#1}}
       V_2 = \begin{array}{r}
               \begin{array}{r}
                 N_c                          \left\{\vphantom{\begin{array}{l}1
                   \\\vdots \\1\\ \end{array}}\right. \\
                 N_c + 1                      \left\{\vphantom{\begin{array}{l}0
                   \\1\\\vdots \\1\\ \end{array}}\right. \\
                 2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}0
                   \\\vdots \\0\\ \end{array}}\right.\\
               \end{array}
               \begin{pmatrix}
               1 \\
               \vdots \\
               1 \\
               0 \\
               1 \\
               \vdots \\
               1 \\
               0 \\
               \vdots \\
               0
               \end{pmatrix}
             \end{array}

    by applying the operation

    .. math::

       V_{H_2}\vert x \rangle \vert 0 \rangle \mapsto
       \left\{
       \begin{array}{lr}
       \vert x \rangle \vert 1 \rangle & \text{if } (x < 2N_c) \land (x \neq N_c) \\
       \vert x \rangle \vert 0 \rangle & \text{else} \\
       \end{array}.
       \right.


    The signs oracle encodes a unconditionally positive sign (as :math:`H_2` has no
    negative entry) by applying the operation

    .. math::

       S_{H_2}\vert x \rangle \vert 0 \rangle \mapsto \vert x \rangle \vert 0 \rangle

    In other words, :math:`S_{H_2}` is the identity operation.

    .. note::

       The :math:`+` sign is encoded as :math:`\vert 0 \rangle`.

    :param n: Size of the :math:`\ket{x}` and :math:`\ket{m}` quantum
        registers.
    :param discretisation_points_number: Number of discretisation points used to
        discretise the wave equation.
    :return: a QRoutine implementing the oracle of the second matrix of the
        decomposition.
    """
    rout = QRoutine()
    x = rout.new_wires(n)
    m = rout.new_wires(n)
    v = rout.new_wires(1)
    s = rout.new_wires(1)  # Unused because it should not be affected by the oracle
    a = rout.new_wires(4)
    rout.set_ancillae(a)

    # Declaring an AbstractGate that will compare the value stored in a quantum
    # register to a constant value.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to compare.
    # comparator = AbstractGate(
    #     "compare_const", [int, int], arity=lambda q_size, rhs: 2 * q_size
    # )
    comparator = AbstractGate(
        "compare_const", [int, int], arity=lambda q_size, rhs: q_size + 1
    )

    # Declaring an AbstractGate that will add a constant value to the value stored in a
    # quantum register.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to add.
    adder = AbstractGate("add_const", [int, int], arity=lambda q_size, rhs: q_size)

    # Declaring an AbstractGate that will subtract a constant value to the value stored
    # in a quantum register.
    # The first int is the size of the input quantum state (number of qubits) and the
    # second int is the constant number to subtract.
    subtractor = AbstractGate("sub_const", [int, int], arity=lambda q_size, rhs: q_size)

    #########################
    # 1. Permutation oracle #
    #########################
    # Copy the value stored in |x> to |m>.
    for i in range(n):
        rout.apply(CNOT, x[i], m[i])

    # |a[0]>  :=  (|x> < discretisation_points_number-2)
    rout.protected_apply(comparator(n, discretisation_points_number - 1), x, a[0])

    # Perform the addition when (|x> < discretisation_points_number-2)
    rout.apply(adder(n, discretisation_points_number - 1).ctrl(), a[0], m)
    # Else, perform a subtraction as the inverse of an addition.
    rout.protected_apply(X, a[0])
    rout.apply(subtractor(n, discretisation_points_number - 1).ctrl(), a[0], m)
    rout.protected_apply(X, a[0])

    # /!\ DO NOT uncompute the comparison as we can re-use it for the weight oracle
    # /!\ part. The compare operation will be uncomputed after the weight oracle.
    # rout.protected_apply(
    #     comparator(n, discretisation_points_number - 1).dag(), x, a[0], m[:-1]
    # )

    #####################
    # 2. Weights oracle #
    #####################
    # Initialise |w> with ones on the 2*(N_d-2)+1 = 2*N_c+1 first points
    # ------------------------------------------------------------------
    rout.apply(comparator(n, 2 * (discretisation_points_number - 2) + 1), x, v)

    # Set |w> to |0> when x == N_c == N_d-2
    # -------------------------------------
    # |a[0]>  :=  (|x> < discretisation_points_number - 2 + 1)
    # /!\ the line below has already been computed by the permutation oracle /!\
    # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    # rout.protected_apply(
    #     arithmetic_compare(n, discretisation_points_number - 1), x, a[0], m[:-1]
    # )
    # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    # |a[1]>  :=  ! (|x>  <   discretisation_points_number - 2)
    #         :=    (|x>  >=  discretisation_points_number - 2)
    rout.protected_apply(comparator(n, discretisation_points_number - 2), x, a[1])
    rout.protected_apply(X, a[1])

    # |v>  :=  ( |a[0]>  &&  |a[1]> )
    rout.apply(CCNOT, a[0], a[1], v)

    # Uncompute
    rout.protected_apply(X.dag(), a[1])
    rout.protected_apply(comparator(n, discretisation_points_number - 2).dag(), x, a[1])
    # /!\ Do not forget the uncomputation below /!\
    rout.protected_apply(comparator(n, discretisation_points_number - 1).dag(), x, a[0])

    ###################
    # 3. Signs oracle #
    ###################
    # We have nothing to do here, all the entries of this matrix are positive so we
    # should leave |s> is the |0> state.

    return rout


# def _get_oracle_dirichlet1_1c_wave_equation_ncc(
#     n: int, discretisation_points_number: int
# ) -> NamedRoutine:
#     r"""Return the oracle for the first matrix in the decomposition.
#
#     The returned routine needs :math:`2n + 8` qubits organised as follow: ::
#
#         |      |x>      |      |m>      |       |v>        |  |s>  |        |a>        |
#         |   .   .   .   |   .   .   .   |  .   .   .   .   |   .   |   .   .   .   .   |
#         |   0      n-1  |   n     2*n-1 | 2*n        2*n+1 | 2*n+2 | 2*n+3       2*n+8 |
#
#     The quantum registers :math:`\ket{m}`, :math:`\ket{v}` and :math:`\ket{s}` should
#     all be initialised to :math:`\ket{0}`.
#
#     This oracle is split in 3 different oracles: an oracle to encode the positions of
#     the non-zero entries (called the permutation oracle), a second oracle to encode
#     the absolute values of the weights of the non-zero entries (called the weights
#     oracle) and a last oracle encoding the signs of the non-zero entries (the sign
#     oracle).
#
#
#     The permutation oracle encodes the permutation
#
#     .. math::
#        \newcommand{\mbf}[1]{\mathbf{#1}}
#        \newcommand{\HmatrixOnePerm}{%
#            \begin{matrix}%
#                0      & \cdots & \cdots & 0      & 1      & 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                \vdots &        &        & \vdots & 0      & \ddots & \ddots &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & 0      & 0      & \cdots & 0      & 1      & 0      & \vdots &        &        &        &        & \vdots \\%
#                1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \vdots &        &        &        &        & \vdots \\%
#                0      & \ddots & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        & \ddots & 1      & 0      &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                0      & \cdots & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                \vdots &        &        &        &        & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        &        & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        &        &        & \ddots & \mbf{1}& 0      &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        &        &        &        & 0      & \mbf{1}& \ddots &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & 0      & \mbf{1}& 0      & 0      & \cdots & 0      \\%
#            \end{matrix}%
#        }%
#        H_1^M = \begin{array}{r}
#                \begin{array}{cccp{5pt}}
#                  \overbrace{\hphantom{\BTmatrix}}^{N_c}&
#                  \overbrace{\hphantom{\Bmatrix}}^{N_c+1}&
#                  \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q-\left(2N_c+1\right)} &
#                \end{array}\\[0pt]
#                \begin{array}{r}
#                  N_c                          \left\{\vphantom{\Bmatrix}\right.  \\[5pt]
#                  N_c + 1                      \left\{\vphantom{\BTmatrix}\right. \\[5pt]
#                  2^q - \left( 2N_c+1\right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
#                \end{array}
#                \left( \HmatrixOnePerm \right)
#              \end{array}
#
#     by applying the operation
#
#     .. math::
#
#        M_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
#        \left\{
#        \begin{array}{lr}
#        \vert x \rangle \vert x + N_c \rangle & \text{if } x < N_c \\
#        \vert x \rangle \vert x - N_c \rangle & \text{else} \\
#        \end{array}.
#        \right.
#
#
#
#     The weights oracle encodes the weights given by the following vector:
#
#     .. math::
#
#        \newcommand{\mbf}[1]{\mathbf{#1}}
#        V_1 = \begin{array}{r}
#                \begin{array}{r}
#                  N_c                          \left\{\vphantom{\begin{array}{l}1
#                    \\\vdots \\1\\ \end{array}}\right. \\
#                  N_c + 1                      \left\{\vphantom{\begin{array}{l}1
#                    \\\vdots \\1\\ \end{array}}\right. \\
#                  2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}0
#                    \\\vdots \\0\\ \end{array}}\right.\\
#                \end{array}
#                \begin{pmatrix}
#                1 \\
#                \vdots \\
#                1 \\
#                1 \\
#                \vdots \\
#                1 \\
#                0 \\
#                \vdots \\
#                0
#                \end{pmatrix}
#              \end{array}
#
#     by applying the operation
#
#     .. math::
#
#        V_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
#        \left\{
#        \begin{array}{lr}
#        \vert x \rangle \vert 1 \rangle & \text{if } x < 2N_c \\
#        \vert x \rangle \vert 0 \rangle & \text{else} \\
#        \end{array}.
#        \right.
#
#
#
#     The signs oracle encodes the signs given by the following vector:
#
#     .. math::
#
#        \newcommand{\mbf}[1]{\mathbf{#1}}
#        S_1 = \begin{array}{r}
#                \begin{array}{r}
#                  N_c                          \left\{\vphantom{\begin{array}{l}+\\-
#                    \\\vdots \\-\\ \end{array}}\right. \\
#                  N_c + 1                      \left\{\vphantom{\begin{array}{l}+\\-
#                    \\\vdots \\-\\ \end{array}}\right. \\
#                  2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}+/-
#                    \\\vdots \\+/-\\ \end{array}}\right.\\
#                \end{array}
#                \begin{pmatrix}
#                + \\
#                - \\
#                \vdots \\
#                - \\
#                + \\
#                - \\
#                \vdots \\
#                - \\
#                +/- \\
#                \vdots \\
#                +/-
#                \end{pmatrix}
#              \end{array}
#
#     by applying the operation
#
#     .. math::
#
#        S_{H_1}\vert x \rangle \vert 0 \rangle \mapsto
#        \left\{
#        \begin{array}{lr}
#        \vert x \rangle \vert 0 \rangle & \text{if } (x = 0) \lor (x = N_c)\\
#        \vert x \rangle \vert 1 \rangle & \text{else} \\
#        \end{array}.
#        \right.
#
#     .. note::
#
#        The :math:`+` sign is encoded as :math:`\vert 0 \rangle`, the :math:`-` sign
#        as :math:`\vert 1 \rangle` and :math:`+/-` can be any of :math:`+` or :math:`-`.
#
#
#     :param n: Size of the :math:`\ket{x}` and :math:`\ket{m}` quantum
#         registers.
#     :param discretisation_points_number: Number of discretisation points used to
#         discretise the wave equation.
#     :return: a QRoutine implementing the oracle of the first matrix of the
#         decomposition.
#     """
#     x = list(range(n))
#     m = list(range(n, 2 * n))
#     v = list(range(2 * n, 2 * n + 2))
#     s = 2 * n + 2
#     a = list(range(2 * n + 3, 2 * n + 9))
#
#     rout = NamedRoutine("oracle_dirichlet1_1D_wave_equation_ncc", arity=2 * n + 9)
#
#     #########################
#     # 1. Permutation oracle #
#     #########################
#     # Copy the value stored in |x> to |m>.
#     for i in range(n):
#         rout.apply(CNOT, x[i], m[i])
#
#     # Compare |x> with a generation-time value.
#     rout.protected_apply(compare(n, discretisation_points_number - 2), x, a[1], a[0])
#
#     # Perform the addition when (|x> < discretisation_points_number-2)
#     rout.apply(add_const(n, discretisation_points_number - 2).ctrl(), a[0], m)
#     # Else, perform a subtraction as the inverse of an addition.
#     rout.protected_apply(X, a[0])
#     rout.apply(sub_const(n, discretisation_points_number - 2).ctrl(), a[0], m)
#     rout.protected_apply(X, a[0])
#
#     # /!\ DO NOT uncompute the comparison as we can re-use it for the sign oracle part.
#     # /!\ The compare operation will be uncomputed after the sign oracle.
#     # rout.protected_apply(compare(n, discretisation_points_number - 2).dag(),
#     #                      x, a[1], a[0])
#
#     ###################
#     # 2. Signs oracle #
#     ###################
#     # |a[1]> := (|x>  <  discretisation_points_number-1)
#     rout.protected_apply(compare(n, discretisation_points_number - 1), x, a[4], a[1])
#
#     # |a[0]> := ! (|x>  <  discretisation_points_number-2)
#     # /!\ the line below has already been computed by the permutation oracle /!\
#     # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
#     # rout.protected_apply(compare(n, discretisation_points_number - 2), x, a[4], a[0])
#     # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
#     # Still, we need the negated value.
#     rout.protected_apply(X, a[0])
#
#     # |a[2]> := |a[0]> && |a[1]>  <==>  (|x> == discretisation_points_number-2)
#     rout.protected_apply(CCNOT, a[0], a[1], a[2])
#     # |a[3]> := (|x>  <  1)  <==>  (|x> == 0)
#     rout.protected_apply(compare(n, 1), x, a[4], a[3])
#
#     # |s>    := ! (|a[2]> || |a[3]>)
#     rout.apply(logical_or, a[2], a[3], s)
#     rout.apply(X, s)
#
#     # Uncompute everything except |s>
#     rout.protected_apply(compare(n, 1).dag(), x, a[4], a[3])
#     rout.protected_apply(CCNOT.dag(), a[0], a[1], a[2])
#     rout.protected_apply(X.dag(), a[0])
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 1).dag(), x, a[4], a[1]
#     )
#     # /!\ Do not forget the uncomputation below /!\
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 2).dag(), x, a[4], a[0]
#     )
#
#     #####################
#     # 3. Weights oracle #
#     #####################
#
#     # Let N := discretisation_points_number
#     # If (|x>  <  (N-2) // 2) OR
#     # (|x>  <  (N - 2 + (N-2) // 2) AND ! (|x>  <  N-2))
#     # Then |v> := 1
#
#     # |a[0]> := (|x>  <  (discretisation_points_number - 2) // 2)
#     rout.protected_apply(
#         compare(n, (discretisation_points_number - 2) // 2), x, a[1], a[0]
#     )
#     # |a[1]> := ! (|x>  <  discretisation_points_number - 2)
#     rout.protected_apply(compare(n, discretisation_points_number - 2), x, a[2], a[1])
#     rout.protected_apply(X, a[1])
#     # |a[2]> := (|x>  <  (discretisation_points_number - 2
#     #                     + (discretisation_points_number - 2) // 2))
#     rout.protected_apply(
#         compare(
#             n,
#             discretisation_points_number - 2 + (discretisation_points_number - 2) // 2,
#         ),
#         x,
#         a[3],
#         a[2],
#     )
#     # |a[3]> := |a[1]> AND |a[2]>
#     rout.protected_apply(CCNOT, a[1], a[2], a[3])
#     # |v[1]> := |a[3]> OR |a[0]> (big-endian for |v> = |01> if condition satisfied).
#     rout.apply(logical_or, a[0], a[3], v[1])
#     # Uncompute a[3]:
#     rout.protected_apply(CCNOT.dag(), a[1], a[2], a[3])
#     # The values in the qubits a[0], a[1] and a[2] will be reused.
#
#     # Let N := discretisation_points_number
#     # If ((|x>  >=  (N - 2) // 2) AND (|x>  <  N - 2)) OR
#     # (|x>  >=  (N - 2 + N // 2) AND (|x>  <  2*N-4))
#     # Then |v> := 2
#
#     # |a[0]> := ! (|x>  <  (discretisation_points_number - 2) // 2)
#     rout.protected_apply(X, a[0])
#     # |a[1]> := (|x>  <  discretisation_points_number - 2)
#     rout.protected_apply(X, a[1])
#     # |a[2]> := (|x>  >=  (discretisation_points_number - 2
#     #                      + (discretisation_points_number - 2) // 2))
#     rout.protected_apply(X, a[2])
#     # |a[3]> := (|x>  <  2 * discretisation_points_number - 4)
#     rout.protected_apply(
#         compare(n, 2 * discretisation_points_number - 4), x, a[4], a[3]
#     )
#     # |a[4]> := (|a[0]>  AND  |a[1]>)
#     rout.protected_apply(CCNOT, a[0], a[1], a[4])
#     # |a[5]> := (|a[2]>  AND  |a[3]>)
#     rout.protected_apply(CCNOT, a[2], a[3], a[5])
#     # |v[0]> := |a[4]> OR |a[5]> (big-endian for |v> = |10> if condition satisfied).
#     rout.apply(logical_or, a[4], a[5], v[0])
#
#     # Uncomputation of everything
#     rout.protected_apply(CCNOT.dag(), a[2], a[3], a[5])
#     rout.protected_apply(CCNOT.dag(), a[0], a[1], a[4])
#     rout.protected_apply(
#         compare(n, 2 * discretisation_points_number - 4).dag(), x, a[4], a[3]
#     )
#     rout.protected_apply(X.dag(), a[2])
#     rout.protected_apply(X.dag(), a[1])
#     rout.protected_apply(X.dag(), a[0])
#     rout.protected_apply(
#         compare(
#             n,
#             discretisation_points_number - 2 + (discretisation_points_number - 2) // 2,
#         ).dag(),
#         x,
#         a[3],
#         a[2],
#     )
#     rout.protected_apply(X.dag(), a[1])
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 2).dag(), x, a[2], a[1]
#     )
#     rout.protected_apply(
#         compare(n, (discretisation_points_number - 2) // 2).dag(), x, a[1], a[0]
#     )
#
#     return rout
#
#
# def _get_oracle_dirichlet2_1d_wave_equation_ncc(
#     n: int, discretisation_points_number: int
# ) -> NamedRoutine:
#     r"""Return the oracle for the second matrix in the decomposition.
#
#     The returned routine needs :math:`2n + 5` qubits organised as follow: ::
#
#         |      |x>      |      |m>      |       |v>        |  |s>  |        |a>        |
#         |   .   .   .   |   .   .   .   |  .   .   .   .   |   .   |   .   .   .   .   |
#         |   0      n-1  |   n     2*n-1 | 2*n        2*n+1 | 2*n+2 | 2*n+3       2*n+7 |
#
#     The quantum registers :math:`\ket{m}`, :math:`\ket{v}` and :math:`\ket{s}` should
#     all be initialised to :math:`\ket{0}`.
#
#     This oracle is split in 3 different oracles: an oracle to encode the positions of
#     the non-zero entries (called the permutation oracle), a second oracle to encode
#     the absolute values of the weights of the non-zero entries (called the weights
#     oracle) and a last oracle encoding the signs of the non-zero entries (the sign
#     oracle).
#
#
#     The permutation oracle encodes the permutation
#
#     .. math::
#
#        \newcommand{\mbf}[1]{\mathbf{#1}}
#        \newcommand{\HmatrixTwo}{%
#            \begin{matrix}%
#                0      & \cdots & \cdots & 0      & 0      & 1      & 0      & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                \vdots &        &        & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        & \vdots & \vdots &        & \ddots & \ddots & 0      & \vdots &        &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & 0      & 1      & 0      &        &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & 0      & \mbf{1}&        &        &        &        & \vdots \\%
#                1      & \ddots &        & \vdots & \vdots &        &        &        & \vdots & 0      &        &        &        &        & \vdots \\%
#                0      & \ddots & \ddots & \vdots & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots & \ddots & \ddots & 0      & \vdots &        &        &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                0      & \cdots & 0      & 1      & 0      & \cdots & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                0      & \cdots & \cdots & 0      & \mbf{1}& 0      & \cdots & \cdots & 0      & 0      & \cdots & \cdots & \cdots & \cdots & 0      \\%
#                \vdots &        &        &        & \ddots & \ddots & \ddots &        & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        & \ddots & \ddots & \ddots & \vdots & \vdots &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        &        & \ddots & \ddots & 0      & 0      &        &        &        &        & \vdots \\%
#                \vdots &        &        &        &        &        &        & \ddots & \mbf{1}& 0      & \ddots &        &        &        & \vdots \\%
#                0      & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0      & \mbf{1}& 0      & 0      & \cdots & \cdots & 0      \\%
#            \end{matrix}%
#        }%
#         H_2^M = \begin{array}{r}
#                \begin{array}{cccp{5pt}}
#                  \overbrace{\hphantom{\BTmatrix}}^{N_c} &
#                  \overbrace{\hphantom{\Bmatrix}}^{N_c+1}&
#                  \overbrace{\hphantom{\RightTopZeroMatrix}}^{2^q - \left(2N_c+1\right)}&
#                \end{array}\\[0pt]
#                \begin{array}{r}
#                  N_c                          \left\{\vphantom{\Bmatrix}\right.  \\[5pt]
#                  N_c + 1                      \left\{\vphantom{\BTmatrix}\right. \\[5pt]
#                  2^q - \left(2N_c+1\right)\left\{\vphantom{\BotLeftZeroMatrix}\right.\\
#                \end{array}
#                \left( \HmatrixTwo \right)
#              \end{array}
#
#     by applying the operation
#
#     .. math::
#
#        M_{H_2}\vert x \rangle \vert 0 \rangle \mapsto
#        \left\{
#        \begin{array}{lr}
#        \vert x \rangle \vert x + (N_c + 1) \rangle & \text{if } x < (N_c + 1) \\
#        \vert x \rangle \vert x - (N_c + 1) \rangle & \text{else} \\
#        \end{array}.
#        \right.
#
#
#
#     The weights oracle encodes the weights given by the following vector:
#
#     .. math::
#
#        \newcommand{\mbf}[1]{\mathbf{#1}}
#        V_2 = \begin{array}{r}
#                \begin{array}{r}
#                  N_c                          \left\{\vphantom{\begin{array}{l}1
#                    \\\vdots \\1\\ \end{array}}\right. \\
#                  N_c + 1                      \left\{\vphantom{\begin{array}{l}0
#                    \\1\\\vdots \\1\\ \end{array}}\right. \\
#                  2^q - \left( 2N_c + 1 \right)\left\{\vphantom{\begin{array}{l}0
#                    \\\vdots \\0\\ \end{array}}\right.\\
#                \end{array}
#                \begin{pmatrix}
#                1 \\
#                \vdots \\
#                1 \\
#                0 \\
#                1 \\
#                \vdots \\
#                1 \\
#                0 \\
#                \vdots \\
#                0
#                \end{pmatrix}
#              \end{array}
#
#     by applying the operation
#
#     .. math::
#
#        V_{H_2}\vert x \rangle \vert 0 \rangle \mapsto
#        \left\{
#        \begin{array}{lr}
#        \vert x \rangle \vert 1 \rangle & \text{if } (x < 2N_c) \land (x \neq N_c) \\
#        \vert x \rangle \vert 0 \rangle & \text{else} \\
#        \end{array}.
#        \right.
#
#
#     The signs oracle encodes a unconditionally positive sign (as :math:`H_2` has no
#     negative entry) by applying the operation
#
#     .. math::
#
#        S_{H_2}\vert x \rangle \vert 0 \rangle \mapsto \vert x \rangle \vert 0 \rangle
#
#     In other words, :math:`S_{H_2}` is the identity operation.
#
#     .. note::
#
#        The :math:`+` sign is encoded as :math:`\vert 0 \rangle`.
#
#     :param n: Size of the :math:`\ket{x}` and :math:`\ket{m}` quantum
#         registers.
#     :param discretisation_points_number: Number of discretisation points used to
#         discretise the wave equation.
#     :return: a QRoutine implementing the oracle of the second matrix of the
#         decomposition.
#     """
#     x = list(range(n))
#     m = list(range(n, 2 * n))
#     v = list(range(2 * n, 2 * n + 2))
#     # The sign qubit is not used by this oracle. In order to save resources, it will
#     # be used as an ancilla qubit (and so returned in the state |0>).
#     # s = 2 * n + 2
#     a = list(range(2 * n + 2, 2 * n + 8))
#
#     rout = NamedRoutine("oracle_dirichlet2_1D_wave_equation_ncc", arity=2 * n + 8)
#
#     #########################
#     # 1. Permutation oracle #
#     #########################
#     # Copy the value stored in |x> to |m>.
#     for i in range(n):
#         rout.apply(CNOT, x[i], m[i])
#
#     # |a[0]>  :=  (|x> < discretisation_points_number-2)
#     rout.protected_apply(compare(n, discretisation_points_number - 1), x, a[2], a[0])
#
#     # Perform the addition when (|x> < discretisation_points_number-2)
#     rout.apply(add_const(n, discretisation_points_number - 1).ctrl(), a[0], m)
#     # Else, perform a subtraction as the inverse of an addition.
#     rout.protected_apply(X, a[0])
#     rout.apply(sub_const(n, discretisation_points_number - 1).ctrl(), a[0], m)
#     rout.protected_apply(X, a[0])
#     # Uncompute
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 1).dag(), x, a[2], a[0]
#     )
#
#     #####################
#     # 2. Weights oracle #
#     #####################
#
#     # Let N := discretisation_points_number
#     # If (|x>  <  (N-2) // 2) OR
#     # (|x>  <  (N - 1 + (N-2) // 2) AND ! (|x>  <  N-1))
#     # Then |v> := 1
#
#     # |a[0]> := (|x>  <  (discretisation_points_number - 2) // 2)
#     rout.protected_apply(
#         compare(n, (discretisation_points_number - 2) // 2), x, a[5], a[0]
#     )
#     # |a[1]> := ! (|x>  <  discretisation_points_number - 1)
#     rout.protected_apply(compare(n, discretisation_points_number - 1), x, a[5], a[1])
#     rout.protected_apply(X, a[1])
#     # |a[2]> := (|x>  <  (discretisation_points_number - 1
#     #                     + (discretisation_points_number - 2) // 2))
#     rout.protected_apply(
#         compare(
#             n,
#             discretisation_points_number - 1 + (discretisation_points_number - 2) // 2,
#         ),
#         x,
#         a[5],
#         a[2],
#     )
#     # |a[3]> := |a[1]> AND |a[2]>
#     rout.protected_apply(CCNOT, a[1], a[2], a[3])
#     # |v[1]> := |a[3]> OR |a[0]> (big-endian for |v> = |01> if condition satisfied).
#     rout.apply(logical_or, a[0], a[3], v[1])
#
#     # Uncompute a[3] and a[1]:
#     rout.protected_apply(CCNOT.dag(), a[1], a[2], a[3])
#     rout.protected_apply(X.dag(), a[1])
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 1).dag(), x, a[5], a[1]
#     )
#     # The values in the qubits a[0] and a[2] will be reused.
#
#     # Let N := discretisation_points_number
#     # If ((|x>  >=  (N - 2) // 2) AND (|x>  <  N - 2)) OR
#     # (|x>  >  (N - 2 + N // 2) AND (|x>  <  2*N-3))
#     # Then |v> := 2
#
#     # |a[0]> := ! (|x>  <  (discretisation_points_number - 2) // 2)
#     rout.protected_apply(X, a[0])
#     # |a[1]> := (|x>  <  discretisation_points_number - 2)
#     rout.protected_apply(compare(n, discretisation_points_number - 2), x, a[5], a[1])
#     # |a[2]> := (|x>  >  (discretisation_points_number - 2
#     #                     + (discretisation_points_number - 2) // 2))
#     rout.protected_apply(X, a[2])
#     # |a[3]> := (|x>  <  2 * discretisation_points_number - 3)
#     rout.protected_apply(
#         compare(n, 2 * discretisation_points_number - 3), x, a[5], a[3]
#     )
#     # |a[4]> := (|a[0]>  AND  |a[1]>)
#     rout.protected_apply(CCNOT, a[0], a[1], a[4])
#     # |a[5]> := (|a[2]>  AND  |a[3]>)
#     rout.protected_apply(CCNOT, a[2], a[3], a[5])
#     # |v[0]> := |a[4]> OR |a[5]> (big-endian for |v> = |10> if condition satisfied).
#     rout.apply(logical_or, a[4], a[5], v[0])
#
#     # Uncomputation of everything
#     rout.protected_apply(CCNOT.dag(), a[2], a[3], a[5])
#     rout.protected_apply(CCNOT.dag(), a[0], a[1], a[4])
#     rout.protected_apply(
#         compare(n, 2 * discretisation_points_number - 3).dag(), x, a[5], a[3]
#     )
#     rout.protected_apply(X.dag(), a[2])
#     rout.protected_apply(
#         compare(n, discretisation_points_number - 2).dag(), x, a[5], a[1]
#     )
#     rout.protected_apply(X.dag(), a[0])
#     rout.protected_apply(
#         compare(
#             n,
#             discretisation_points_number - 1 + (discretisation_points_number - 2) // 2,
#         ).dag(),
#         x,
#         a[5],
#         a[2],
#     )
#     rout.protected_apply(
#         compare(n, (discretisation_points_number - 2) // 2).dag(), x, a[5], a[0]
#     )
#     ###################
#     # 3. Signs oracle #
#     ###################
#     # We have nothing to do here, all the entries of this matrix are positive so we
#     # should leave |s> is the |0> state.
#
#     return rout
